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Stochastic dynamics of N bistable elements with global time-delayed interactions: 
towards an exact solution of the master equations for finite N 
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We consider a network of N noisy bistable elements with global time-delayed couplings. In a two- 
state description, where elements are represented by Ising spins, the collective dynamics is described 
by an infinite hierarchy of coupled master equations which was solved at the mean-field level in the 
thermodynamic limit. For a finite number of elements, an analytical description was deemed so far 
intractable and numerical studies seemed to be necessary. In this paper we consider the case of two 
interacting elements and show that a partial analytical description of the stationary state is possible 
if the stochastic process is time-symmetric. This requires some relationship between the transition 
rates to be satisfied. 



^^ PACS numbers: 02.50.Ey, 05.40.-a, 05.10.Gg 

o 

(D I. INTRODUCTION 



I '. Systems with time-delayed interactions have been a subject of extensive studies in recent years due to their relevance 

C^ to a wide range of phenomena occurring in physics, biology, ecology, economics, and other sciences. In most situations, 
c/j the effect of random noise due to the environmental fluctuations cannot be ignored and observables are better described 
^ as stochastic variables. As is well known, this may have a major impact on the dynamical behavior, in particular 

J5 when noise combines with nonlinearity, which leads to many remarkable effects such as dynamical transitions [T], 
JI^ synchronization [2], stochastic resonance^, coherence resonance[l], etc... The addition of time delay increases the 

' , dimensionality and hence the complexity of these systems, inducing new phenomena such as multi-stability[S] and 

^ oscillatory behavior [MSI . 

Q By definition, stochastic time-delay systems are non-Markovian, which seriously complicates the analytical treat- 

^ ment as the standard tools for ordinary stochastic differential equations are not directly applicable. Although one 

'~~' can extend the Fokker-Planck description to stochastic delay-differential equations [3 HH]: exact solutions are rare 

I (essentially limited to the linear case[9j[TTl[T2]), and in order to calculate probability densities and time correlation 

^ functions one has to resort to approximate treatments {e.g. small-delay expansion [Ul I13j. perturbation theory P^) 

C^ and more generally to numerical simulations. This can be traced back to the fact that stochastic delay systems can 

lO be viewed as systems with an infinite number of degrees of freedom. 

'^ The situation is even more complicated when one considers several interacting units with time-delayed couplings, 

^^ a situation that usually occurs in a biological context {e.g. in neural or genetic regulatory networks) and can also be 

f — realized with laser networks (see e.g. [TSl [H] for recent references). From a theoretical point of view, the canonical 

^^ example is a globally coupled network of stochastically driven bistable elements, a system that has been extensively 

^> studied in the absence of delay [Sj [T71 [18]. As is often done with bistable elements, one may replace the original 

\ I continuous system by a two-state model with suitably chosen transition rates and replace the stochastic differential 

>• equations by master equations for the occupation probabilities. Then one has to cope with an infinite hierarchy of 

coupled equations which at first sight cannot be closed. So far, a full analytical description of the dynamics is only 

, . avail able when there is a single element [191 I20j or an infinite number of elements [21]. In this thermodynamic limit, 

^ which may be justified in actual situations {e.g. a multicellular system[8]), one can derive a deterministic equation 

of motion for the mean-field variable i.e. the ensemble average of the state variable. The solution then exhibits 

phase transitions to nontrivial stationary states or delay-dependent oscillations via Hopf bifurcations. In principle, 

corrections to the mean-field behavior can be obtained within an expansion in the inverse system size (see e.g. j22|). 

The aim of the present work is to show that a solution of the delay master equations can also be obtained for a 

finite number of interacting bistable elements, at least in the stationary state and for a restricted time interval. In 

the following, for simplicity, we only treat the case of two coupled elements but the demonstration can be extended to 

several units at the price of increasing analytical complexity. Interestingly, the demonstration used the time-symmetry 

of the delay master equations, an issue which does not seem to have been discussed in the existing literature. 

The paper is organized as follows. In the next Section we present the model and derive the master equations for the 



probabilities. In Section III we solve these equations and compute the time correlation functions under the condition 
that a certain relation between the transition rates is satisfied. Analytical results are then compared to the results 
of numerical simulations. Concluding remarks are given in Section IV. Appendix A is devoted to the analysis of 
time-symmetry and Appendix B details the solution of the master equations. 

II. MODEL AND MASTER EQUATIONS 

As Huber and Tsimring[5T] we consider an ensemble of N identical bistable elements, each of them characterized 
by the variable Xi{t) and obeying the coupled Langevin equations 

^^^^^ = _^KM + ,x{t - r) + ^/2DUt) , I = 1...^ (1) 

dxi 

where V{x) = —x'^/2 + a;*/4 is a generic symmetric double-well potential and X{t) is the global 'mean-field' 

i 

Here r is the time delay, e is the strength of the feedback coupling, and D is the variance of the Gaussian fluctuations, 
which are (5-correlated and mutually independent < £,i{t)£,j{t') >— 6{t — t')5ij. For each i, this set of equations thus 
describes the overdamped motion of a particle evolving in an effective r-dependent double- well potential Ur{x) — 
V{x) — eX{t — t)x. Note that each element is identically coupled to all units at time i — r, including itself (the case of 
a chain with unidirectional coupling, i.e. Xi{t) only coupled to Xi-iit — r), is much simpler and has been considered 
inRef.[13]). 

In the following we shall be interested in the stationary state that is reached in the large-time limit. Since the 
number of elements is strictly finite, one expects this state to be unique for all couplings with an average value 
< X{t) >~ {i.e. there is no phase transition). 

As in Refs.fini [H] we consider the case of small noise and small coupling where one can neglect the intrawell 
fiuctuations and replace the continuous dynamical variables Xi (t) by the two-state variables Si (t) that take the values 
±1. The switching rates associated to the instantaneous potential Ur{x) can be calculated from Kramers formula [24] 
7k = {2'k)~^ \JU!^{x±)U'^{xo) exp{— AUr / D) where x± and xq are the positions of the minima and the maximum of 
the potential, respectively, and AUr = Ut{xo) — Ut{x±) is the barrier that an element has to overcome to jump from 
one stable state to the other. For small e, the two minima of the potential are located at x± = ±1 + eX/2, which 
yields m] 

v/2 ± SeXjt ~ t) l±4eA(t-r) 
7x = Y^ -P( 4^ )' (3) 

so that 7if can take different values depending on the sign of Si at time t and the state of the system at time t ~ t 
(in the two-state description X{£) ~ (1/A^) X^j Si{t)). When A^ = 1, which is the simplest case studied by Tsimring 
and Pikovskv|19j. there are only two possible values 



V2±3e . l±4e. 

71,2 = ^^exp(-^^) (4) 

corresponding to s{t)s{t — t) = 1 and s{t)s{t — t) = —1, respectively. In the opposite limit studied in Ref. [21, where 
the number of units is very large, stochastic fluctuations can be neglected and the collective variable X{t) approaches 
the average value < s{t) >— J2 s=±i ^Pi^ ' ^) > where p(±l,i) are the occupation probabilities of the states s = ±1. 
One can then derive a closed equation of motion for the mean field X{t). 

In the present work we focus on the case N — 2 and from Eq. (Isl) we must take into account three different switching 
rates 



70 = ^^^^exp(-i^) for |A(t - t)| = 1 and s,(i)A(i - r) = 1 

ZTT 41^ 

/2 1 

71 = |^°^P(-i;^) forA(i-T)=0 

72 = ^^^^exp(-i^) for |A(t - r)| = 1 and s,(t)A(i - r) = -1 (5) 



with X{t) — {1 /2)[si{t) + S2{t)]. These different rates can be put together in a single expression defining the switching 
rate from Si{t) to —Si(t) depending on the state of the spins si and S2 at t — r, 

T{s, -^ -~s,;s) = 71 + 1(70 - J2)s^{t)[sl{t ~ t) + S2{t - r)] + ^(70 + 72 - 27i)[l + si{t - T)s2{t ~ r)] (6) 

where s denotes the two-component vector {si,S2} and it is imphcit in the notation T(si — > — Si;s) that Si is taken 
at time t and s at time t — t. 

The description of the dynamics of the system is then encoded in 4 coupled master equations for the probabilities 
p{s, t) which read 

p{{si,S2}.,t) =^ \t{-si -^ si;s')p(s',t-T;{-si,S2},t)+r(-S2 -^ S2;s')p(s',t-T;{si,-S2},t) 

- [r(si -^ -si; s') + r(s2 -^ -S2; s')Ms', t - r; {si, S2}, i)] (7) 

where p{s',t — r;s,t) is the joint probability that the system is in state s' dX t — t and s at t. The equations of 
motion are thus not closed at the level of one-time probabilities and one also need to consider the equations of 
motion for p{s',t — T;s,t), p{s",t — 2T;s',t — r;s,i), and so on. As already pointed out, this hierarchical structure 
merely reveals that a time-delayed system is a system with an infinite number of degrees of freedom. For TV = 1, 
the master equations for the one-time probabilities p{±l,t) are closed[TH] because one can use the exact relations 
p(l, t - r; ±1, i) +p{-l, t-T;±l,t) = p{±l, t) and p{±l,t - r; 1, i) +p{±l, i - r; -1, i) = p{±l, t-r) to eliminate the 
two-time probabilities. For similar reasons, a closed chain with unidirectional couplings can also be solved at the level 
of the one-time probabilities |23j . On the other hand, for the system described by the coupled Langevin equations (II]), 
one can easily check that there are not enough such exact relations to close the hierarchy for a generic value of N, 
even when taking into account the additional symmetry relations obtained by exchanging the spins 1 and 2 and the 
signs -1-1 and —1 (see below). 

III. STATIONARY SOLUTION OF THE MASTER EQUATIONS FOR iV = 2 

We now turn our attention to the stationary solution of the N = 2 model which satisfies pst{s,t) = 0. We are 
interested in calculating Pst(s) and the self and cross time-correlation functions tpsit) and 4'cit) defined by 

Mt) =< si{t)s,{0) >= ^ sis? p,t{s, t\s°) p,t{s°) 
V'c(t)=<si(i)s2(0)>=^sis§p,t(s,t|s0)p,t(s°) (8) 

where we have introduced the conditional probabilities Psf(s,i|s°) (hereafter it is implicit that s'~' denotes the state at 
t = 0). From Eqs. (I?]), we readily get 

p(s, t\s") ^ Yl [^(-*i ^ 5i; s'Ms', t - r; {-si, S2}, t|s°) + T(-S2 ^ S2;s')p{s\ t - r; {si, -S2}, t|s°) 

- [r(si ^ -si; s') + T{s2 ^ -S2; s')]p(s', t - r; {si, S2}, i|s°)] (9) 

which, unsurprisingly, indicates that the calculation of ips{i) and ^c(i) requires the knowledge of the three-time 
conditional probabilities Pst(s', i — t; s, t|s°) which in turn depend on four-time functions, etc... 

Hence, at this stage, it seems that the problem cannot be solved exactly. However, remarkably, an analytical 
solution does exist if the switching rates satisfy the relation 

7072 = 7? (10) 

which may also be viewed as a natural consequence of Kramers' equations (Is]) if one expands the product 7072 up to 
order e. Indeed, as shown in Appendix A, the stochastic process is then statistically time reversible in the stationary 
state. In other words, any sequence of states S has the same probability as its time reverse 5*: 

Pst{S) = p,t{S) . (11) 



This readily implies that 



and thus 



Pst{s', t - r; s°, 0; s, t) = pst{s, -t; s°, 0; s', t - t) 



Pstis', t-T;s, t\s°) = pst{s, -t; s', r - t|s°) 



(12) 



(13) 



(In Appendix A we show an example where Eq.( 10 ) is violated and thus Eq. ( [T3| ) is not satisfied). From this equation 
we can derive a closed equation for the key quantity Pst(s', i — r; s, i|s'^) and then compute the stationary occupation 
probabilities and the time correlation functions. We stress, however, that the solution for the time correlation functions 
is only valid in the interval < i < t (whereas Eq. (Ils]) holds for all times t). In this sense the problem is only 
partially solved. 
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FIG. 1: Sequence of states s, s' and s" for (a) < i < r and (b) t > t. On the right hand side, the direction of time is reversed. 



We first compute the infinitesimal change dp{s' , t — t;s, t\s^) associated to an infinitesimal change dt. Taking into 
account the fact that t appears twice, we obtain 

dp{s', t - r; s, i|s°) = [p(s', t - r; s, i + dt\s°) - p{s' , t - r; s, t|s°)] + [p(s', t + dt - r; s, t|s") - p(s', t - r; s, f |s°)] 

= [p(s', t - r; s, i + dt\s°) - p{s' , t - r; s, t|s°)] + [p(s, -t;s',T-t- dt\s°) - p{s, -t; s', r - t\s°)] 
= dpi + dp2 (14) 



where we have used Eq. (131 to reverse the direction of time in dp2, the second term inside brackets (from now on we 
omit the subscript {st} in the various probability distributions to simplify the notations). When < t < t, the state 
s° of the system at i = is irrelevant for calculating dpi/dt and dp2/dt. Indeed, as illustrated in Fig. 1(a), one has 
t — T<0<iin the first case and — t < < r — i in the second case. Since the switching rate at time t only depends 
on the states at times t and t — t and the switching rate at time t — t only depends on the states at times t — t and 
—t, we simply have (cf. Eq. ([7| without the summations) 

p{s',t-T;s,t\s") =pi +P2 



with 



and 



pi =T{-si -^ si;s')p{s\t - t; {-si, sa}, i|s°) + T{-S2 -^ S2; s')p(s', t - r; {si, -sa}, i|s°) 
- [r(si ^ -si; s') + T{s2 -^ -S2; s')]p(s', t - r; {si, S2}, t|s°) 

P2=- T{-s[ ^ s[;s)p{s, -t; {-s'l, 4}, r - t\s°) - T{-s'^ ^ 4; s)p(s, -t; {4, -4}, r - t|s") 

+ [r(4 ^ -4; s) + r(4 ^ -4; s)]p(s, -t; {4,41, r - t|s") 



(15) 



(16) 



(17) 



(note the change of sign in p2 in relation to the change of sign of dt) . 

On the other hand, when i > r, as illustrated in Fig. 1(b), one has — t < r — t < and a possible switching at time 
T — t is conditioned by the fact that the system is in state s" at i = 0, which lies in the future. Therefore the simple 
probabilistic argument that leads to Eq. (17) is no more valid (for instance, it can be checked in the case A^ = 1 
that the corresponding equation does not yield the correct expression of the correlation function for r < i < 2t as 
computed in Ref.[l9]). Despite our efforts, we have not succeeded so far to find the correct equation of motion of 
Pst(s', i — r; s, i|s°) for t > t. 

To proceed and solve the set of coupled linear differential equations represented by Eqs. (15 171, it is now necessary 
to make explicit the dependence on the 4 configurations {1,1}, {1,-1}, { — 1,1}, { — 1,— l}~that we shall denote 
A, B, C, D, respectively. First, we note that all quantities must be invariant under particle exchange si ^^ S2 and sign 
exchange +1^—1 (this symmetry is unbroken because no phase transition is expected). This readily implies that 
p{B) = p{C) and p{A) = p{D) whence 



piA)+p{B) 



(18) 



Moreover, when calculating the conditional probabilities, the state s" at t = may be chosen to be either A or B. 
Hence there are only 6 distinct functions p(s,i|s°), na.me\y p{A,t\ A), p{B,t\ A), p{D,t\ A), p{A,t\B),p{B,t\B),p{C,t\B) 
and only 3 of them are linearly independent. Indeed, from the equations expressing the conservation of probabilities 



^p(s,t|s") = l 

S 

^p(s,i|sXs«)=p(s), 



(19) 



one can derive the 3 relations 



p{A)p{B,t\A)=p{B)p{A,t\B) 
2p{B,t\A) +p{A,t\A) +p{D,t\A) = 1 
2p{A,t\B) +p{B,t\B)+p{C,t\B) = 1 . 



(20) 



There are also only 2 x 10 distinct functions among the 64 conditional probabilities Pst{s' , t — r; s, i|s°) (which moreover 
are not all linearly independent because of conservation of probabilities) . Inserting the expression mf of the rates in 
Eqs. ( 15p7l, we thus obtain two sets of 10 coupled linear differential equations which may be written as 



PA{t)=TAPA{t) 
PB(t)=rBPB(t) 



(21) 



with 



PA{t) = 



fp{B,t 
piA,t- 
p{D,t 
piCt^ 
p{B,t 
p{A,t- 
p{D,t 
p{B,t^ 
p{A,t- 

\p{D,t 



T-B,t\A)\ 
T;B,t\A) 



B,t\A) 
B,t\A) 
A,t\A) 
A,t\A) 
A,t\A) 
D,t\A) 
T;D,t\A) 
r;D,t\A)J 



, PB{t) = 



fp{A,t~T;A,t\B)\ 
p{B,t-T;A,t\B) 
p{C,t-T;A,t\B) 
p{D,t-T;A,t\B) 
p{A,t~T;B,t\B) 
p{B,t-T-B,t\B) 
p{C,t-T;B,t\B) 
p{A,t-T;C,t\B) 
p(B,t-T;C,t\B) 

\p{C,t^T;C,t\B)J 



(22) 



and Ta, ^b are 10 x 10 matrices whose expressions are given in Appendix B. The solution of these equations is 

PAit) = e^'^'pAiO) 
PBit) = e^'^'psiO) , 



(23) 



and the problem reduces to the calculation of the eigenvalues and eigenvectors of Ta and r^, as detailed in Appendix 
B (the two matrices have the same spectrum when Eq. (10) is satisfied). It only remains to determine the initial 



conditions at i = 0. One can see from Eqs. (22) that only 3 components of the vectors Pa{0) and Pb(0) are nonzero, 



and 



piA, -r; A, 0\A) = p{A, 0; A, t\A) = p(A t\A) 
p{B, -r; A, 0\A) - p(A, 0; B, t\A) = p{B, t\A) 
p{D, -T- A, Q\A) = p{A, 0; D, t\A) = p(D, t\A) 



p{A, -T- B, 0\B) = p{B, 0; A, t|B) = p(A, t|B) 
p(S,-r;B,0|S) =p(S,0;B,r|B) =p(B,t|S) 
p(C, -r; B, 0|S) = p(S, 0; C, t\B) - p(C, t\B) 



(24) 



(25) 



where we have used the time-reversal symmetry, Eq. (13). The 6 quantities p(A,t\A),p{B,t\A), ...,p{C,t\B) are still 
unknown but they can be obtained by observing that Eqs.( 24|25 ) also correspond to the values oi p{s',t — T;s,t\s^) 
at t = T. Therefore, they are also given by 



Pa[t) 
Ps(t) 






(26) 



which yields two sets of self-consistency equations whose solution is given in Appendix B. Interestingly, these equations 
have nontrivial solutions only when Eq. ( 10 ) is satisfied. Therefore this relation again appears as a necessary condition 



for the problem under study to be integrable. 

Knowing these quantities we can calculate all the components of PA{t) and ps(i) (these functions are sums of six 
exponential factors -see Appendix B- and the explicit expressions are not given here for the sake of brevity) , then the 
probabilities p(s,t|s'^) = J2s'P(^''^ ~ ''"! 8,^18"), and finally p{A) and the correlation functions 4's{t),'4'c{t)- We have 
verified that the probabilities p{s, t|s°) obtained in this way satisfy the equations of motion ([9|, which shows that the 
whole calculation is indeed consistent. 

The correlation functions are calculated using 



i^s{t)^2p{A)[p{A,t\A) 
Mt)^MAMA,t\A) 



- p{D, t\A)] + 2p{B) [p{B, t\B) - p{C, t\B)] 
'P{D,t\A)] - 2p{BMB,t\B)-p{C,t\B)] 



(27) 




FIG. 2: Stationary occupation probabilities p{A) and p{B) = 1/2 — p{A) as a function of the time delay r for D = 0.05, 
e = 0.05 (left panel) and e = —0.05 (right panel). Theoretical results (solid lines) are compared to numerical simulations of the 
stochastic process described by the rates ([5| (circles) and simulations of the original Langevin equations ([l| (crosses). 



To test the validity of our analytical results, we compared them to numerical simulations of the stochastic two-state 
process with switching rates obtained from the Kramers relations (Isl) using D = 0.05 and e = ±0.05. For e > this 
yields 70 = 0.000578,72 = 0.003965, and 71 = 0.001516 (for e < 0, the values of 70 and 72 are interchanged). 71 
was actually adjusted to the value ^7072 = 0.001514 so to exactly satisfy Eq. (10). The numerical simulations were 
carried out using a time-step At = 0.01 and averages were taken over a run of 16 x 10^ steps, discarding the first 10^ 
steps. We also performed simulations of the Langevin dynamics described by Eq. (fl]) using Euler method. 

Fig. [2] shows the dependence of the stationary probabilities p{A) and p{B) — 1/2 — p{A) on the time delay r. As 
can be seen, the theory is in perfect agreement with the simulations of the stochastic process. The agreement with 
the Langevin dynamics is also reasonably good for e > whereas some systematic deviations are observed for e < 0. 
There are indeed more fluctuations in the latter case and the random variables Xi{t) and 12 (i) often take intermediate 
values between —1 and -1-1 making the two-state model less accurate. 

It is worth noting that p{A) goes to a nontrivial limit when r — ?► 00 (given analytically by Eq. (B17)). Indeed, the 
trivial value 1/4 would be obtained by naively assuming that the events at t and i — r in the master equations ([7]) 
can be decoupled when r is much larger than the other characteristic times of the system. The actual limit is larger 
than 1/4 and is the same for e > and e < 0, implying that the two configurations A and D where the two spins are 
in the same state are favored whatever the sign of the feedback (this contrasts with the behavior for small time delays 
where the ratio p{A)/p{B) is larger or smaller than 1 depending on the sign of e: this result can be easily understood 
by referring to the case r = where the process becomes Markovian) . To understand the behavior for larger time 
delays, it is instructive to consider the conditional probabihties p{A,t\A),p{B,t\A), ...,p{C,t\B) whose dependence 
on T is shown in Fig. [3] 
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FIG. 3: Probabilities p(A,ri A), p(B,r|A),p(D,T|yl),p(yl,rlB),p(B,rlB),p(C,rlB) as a function of r for e = 0.05 (left panel) 
and e = —0.05 (right panel). 



We note in particular that p(A, t|A) is larger or smaller than p{D,t\A) depending on the sign of e. This is not 
surprising as the global coupling tX(t — r) increases or decreases the probability for an element to be at time t in 
the same potential well in which the majority of elements were at time t — t depending on whether e is positive or 
negativepTI. Since configurations B and C do not contribute to the mean field X, the force that determines the 
evolution of the state A at time t is (on average) e < X(t — t) >— e[p{A,t — T\A,t) — p{D,t — T\A,t)] which is 
equal to e[p(^,T|^) — p{D,t\A)] in the stationary state. This force is thus positive on average whatever the sign 
of e and the configuration A is stabilized. On the other hand, for configuration B, the average force is equal to 
e[p{A,T\B) — p{D,t\B)] which is zero by symmetry. Although this is a mean-field argument, it qualitatively explains 
why configurations A or D are more probable than B or C for both positive and negative feedbacks, as observed in 
Fig. [2] for large r. 



The behavior of the occupation probabihties as a function of the time delay observed in Figs. [2] and Is] (note also 
the non-monotonic variation oip{B, t\A) in this latter figure) suggests that the dependence on the coupling strength 
at fixed t may be also interesting. This is illustrated in Fig. |4] where the ratio p{B)/p{A) computed from Eq. (B15) 
is plotted as a function of e for r = 500. We see that there is a range of negative values of e where p{B) is slightly 
larger than p{A) with a maximum occurring at e ~ —0.012. 




FIG. 4: The ratio of the stationary occupation probabihties p{B)/p{A) as a function of the coupHng strength e for D = 0.05 
and r = 500. This quantity is computed from the analytical expression of p{A) given in Appendix B. 
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FIG. 5: (Color on line) Self (black) and cross (red) time correlation functions ipsit) and Vc(i) for D — 0.05, r = 2000, 
e = 0.05 (left panel) and e — —0.05 (right panel). For < i < r, the theoretical results (solid lines) are compared to numerical 
simulations of the stochastic process (circles). 



Finally, the self and cross time correlation functions for r = 2000 are shown in Fig. [sj For |e| = 0.05, this value of t 
is slightly larger than the largest characteristic time of the system (to = I/70 ~ 1729 for e > 0). There is again perfect 
agreement between theory and simulations of the stochastic two-state process in the interval [0, r] where the analytical 
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results are available (note that V'c(O) = 4p(A) — 1). For e > both functions are positive with maxima at t « nr 
whereas for e < the peaks at t « nr have alternating signs (of course, the functions go to as t — )■ oo). Moreover, 
the peaks are always delayed with respect to nr, a behavior that was already observed in the case N = IfTO], 

IV. SUMMARY AND CONCLUSION 

We have studied a system of two bistable elements with a global time-delayed coupling using a two-state model 
where the dynamics is described by delay-differential master equations. In general, due to the non-Markovian nature 
of the dynamics, this set of equations is not closed since one-time occupation probabilities depend on two-time 
probabilities, two-time probabilities on three-time probabilities, etc... We have shown, however, that one can close 
this infinite hierarchy and derive analytical expressions for the occupation probabilities and the correlation functions 
in the stationary state provided the stochastic process is time-symmetric. This property only holds when some 
relationship between the different switching rates is satisfied, which is approximately the case when the rates are 
described by Kramers's theory in the limit of very small coupling. It is rather obvious that the explicit demonstration 
presented in this work can be generalized to a larger number of interacting elements, which opens the way for a 
full (though admittedly intricate) analytical description. We stress, however, that our solution for A^ = 2 is still 
incomplete since the analytical expressions of the time correlation functions are only valid for t smaller than the time 
delay r, which precludes the calculation of the power spectrum as was done for N = 1[19, or N — > oo[21 . Obtaining 
an analytical description valid for all times is clearly the most challenging task for future work. It would also be 
interesting to compute the distribution of residence times along the lines of Ref . |20j . Our results for TV = 2 show 
that the main effect of the global coupling when the time delay is not too small is to increase the probability for the 
elements to be in the same potential well, whatever the sign of the feedback coupling. As N increases, one should start 
to observe on a certain time-scale the nontrivial behavior exhibited in the thermodynamic limit [51]. For instance, 
since there exists an ordered phase with a non-zero stationary mean field X for a sufficiently large positive coupling, 
one should see X switching between these non-zero values with a rate decreasing with N |25j . 

Appendix A: Time-symmetry in the stationary state 

In this appendix we show that the stochastic process described by the hopping rates ([5]) is statistically time- 



symmetric in the stationary state when relation ( 10 1 is satisfied. Time reversibility is not at all obvious because 
stochastic processes with delay are non-Markovian by definition. However, the probability of observing a given path 
during a time interval [t, t + t] only depends on the realization of the process during the preceding time interval 
[t — T,t] , a property that we shall use repeatedly in the following (more generally, this allows one to describe delay 
processes in terms of Markov processes at the price of enlarging the space of random variables, see e.g. PUIIT^ ). 

For this purpose we consider a discrete time approximation of the original process by defining a set of equidistant 
observation times nAt with a time-step 

Ai = — (Al) 

where M is some large integer. The discrete time process converges to the original continuous one as M — > oo. The 
random path (st„ = s'^''^ s^^ = s'^^\ ..., St,^i — s^*^)) at times to,ti — to + At, ...,tM — to + t is then denoted S^'^^ and 
the conditional probability of observing the path S'M+i,2Af+i ^^^^^ ^^^ p^^j^ go,M jg denoted by p(s^i+h2M+i^go,M^ 
(as we only study the stationary state, we can set to = without loss of generality). We shall first prove that 

p,,(5°'*0=P.t(^'''"), (A2) 

and then extend the proof to a path of arbitrary duration. To simplify the notation, the subscript st is dropped in 
the following and the time reversal of a path S (for instance 5*^'°) is denoted S. 

1. The case TV = 1 

As a warm-up, we first consider the case A^ = 1. As we shall see, time symmetry is always valid, even without 
taking the continuous limit At — > 0. This contrasts with the case N = 2 that is studied in the next subsection. 



To prove Eq. ( A2 ) we consider sample paths of increasing length such that the path in the last time interval of 
duration r is the same as in the first interval [0, r], as illustrated schematically in Fig. [6] We thus first start with a 
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/^ "■■-.. ..-- ^_ 

S k- 1 S2 

t U T 

Sk S 1 Sk Si 



FIG. 6: Description of the sequences of trajectories of increasing length that are used in the proof of Eq. (A2|. In each 



case, the first and last trajectories of duration r are identical. On the right-hand side, the time-reversal of the sequences is 
considered. 



sample path S ^ Sx ^ S ai duration 3t -f 2Ai that reads 5'"'3*^+2 = {st„ = s^^'^.st^ = s'^^\. 



^Stn 



* 7 ^t J 



g{M+l) 



,(A/-|-2) 



, St, 



,(2J\/+1) 



St2 



s(") St 



,(1) 



■•,St3 



a(M)\ 



(here S (resp. Si) is a 
short-hand notation for the sample path in the intervals [0, r] and [2r -f 2Ai, 3t + 2Ai] (resp. [r -I- Ai, 2t -|- Ai])). We 
then consider the time reverse of this path, i.e. S ^ Si ^ S and prove that 

P{Si\S)P{S\Si) = P(S~i\S)P(S\S~i) . 

The crucial point is that the probabihty of observing s'^'^ only depends on s^*^^) 
time step, and s*^'^*^^^\ the value of the spin M + 1 steps earlier. Hence 

2M+1 



(A3) 
the value of the spin at the preceding 



PiSi\s)= n pi^^'^\^ 



Wu(«-i) M-M-1) 



(A4) 



i=M+l 



where the probabihties p{s^''''\s 



(OU(i-l) M-M^l) 



s^ ') are given by t 


he transi 


tion rates c 


p(+l| + !,+!) =p(-l| - 


-1,-1) = 


= 1 - 71 At 


p(-l| + l,+l)=p(+l|- 


-1,-1) = 


= 71 At 


p(+l| + l,-l)=p(-l|- 


-1,+1) = 


= 1 - 72At 


p(-l| + l,-l)-p(+l|- 


-1,+1) = 


= 72At . 



Eq. (A4) can thus be recast as 



P(^i|5)-(l-pi)™^pr(l-P2)™^P^= 



(A5) 



(A6) 



where pi = 7iAt, p2 = 72At, and mi,m2,ni,n2 count how many times the factors 1 — Pi, 1 — P2,Pi,P2 appear in the 
right-hand side of Eq. ( A4 ) . These numbers are readily obtained by introducing the binary variables 



.s^"') + 1 
2 '' 



Ji^+M+l) 



1 



0<iy<M 



(A7) 
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which yields 



M 

mi = 'Y^lauK^iK + (1 - a^){l - 6,y-i)(l - K)] 

M 

ni = ^[a^fe^_i(l - bi,) + (1 - ai,){l - fe^-i)6i.] 

M 



1712 = ^[a^(l - fc,y-i)(l - 6,.) + (1 - a^)&,y-l&i.] 
1^=0 

M 

n2 = Y.^a,^^ - K^i)K + (1 - a.)6.-i(l - h,)] . (A8) 






with the convention 6_i = om- Similarly, we can also write the left and right-hand sides of Eq. (A3| as {\—pi)^p-^{l — 
P2)^'^P2 find (1— pi)'""^p/(l— p2 )^^P 9^ respectively. The explicit calculation readily shows that fci = fci, ^2 = ^2,^1 = Ii 



and ^2 = ^2, which proves Eq. (A3). 

We then intercalate between S and Si an additional path ^2 of duration t, as shown in Fig [61 and we repeat the 
same calculation. After some algebra we find that 

P{S2\S)P{Si\S2)P{S\Si) = P(Sl\S)P{T2\Sl)P{S\S~2) . (A9) 

It is then straightforward to prove that 

P{Sk\S)P{Sk^i\Sk) X • • • X P{S\Si) = P{S[\S)P{S~2\Si) X ... X P(S\S~k) (AlO) 

for an arbitrary integer fc > 1. Summing over the intermediate states (S*!, 82, ■■■, Sk) we then find that 

P{S,k\S)^P(S,k\S) , (All) 

where P{S, k\S) is the conditional probability to realize the path S after k intermediate paths of duration r, given the 
path S in the first time interval [0,r]. Assuming ergodicity, we can forget the initial condition in the limit /c — >■ cx), 



which yields Eq. (A2|. 

The proof can easily be extended to a path of duration nr (and more generally of any duration). Take for instance 
n — 3 and consider the path Sc — >■ Sk — > Sk^i--- — >■ 6*1 — > 5a — )• 5'f, — > Sc where each individual path is of duration r. 



Then from Eq. (AlO I we have 



PiSk\S,)PiSk-l\Sk) X ... X P{Sa\Si)PiSb\Sa)PiSc\Sb) = PiSb\S,)PiSa\Sb)P{Si\Sa) 



x---xPiSk\Sk-i)PiSc\Sk) . (A12) 

Then summing over all intermediate states (^i, 6*2... 5*^) we obtain 

P{Sa, k{T + At)\Sc)P{Sb\Sa)P{Sc\Sb) = P(Sb\S;)P(Sa\Sb)P(S;, k{T + M)\Sa) (A13) 

where P{Sa,k{T + /S.t)\Sc) (resp. P(5'c, fc(T + Ai)|S'a) ) is the conditional probability to reahze the path Sa (resp. Sc) 
after k intervals of duration t + At given the path Sc (resp. So)- In the limit k — > 00, P{Sa, k{T + At)\Sc) — > p{Sa) 
and P(S^, k{T + At)\S^) -^ p(Sc), whence 

p{Sa)P{Sb\Sa)P{Sc\Sb) = P(S~b\S~c)P(S;\S~b)p(S~c) (A14) 

which may be recast as 

p{Sa)P{Sb,Sc\Sa) = P(Sa,Sb\Sc)p(S,) (A15) 

and finally 

p{Sa,Sb, Sc) = p{Sc, Sb, Sa) (A16) 

where {Sa, Sb, Sc) is a path of duration 3t. 
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2. The case iV = 2 



The proof for the case A^ = 2 is conducted along the same Unes, startmg with the sample path S ^ Si ^ S oi 
duration 3t + 2At, where s is now the two-component vector (si,S2) and the rates are given by Eq. ^. Eq. (A4) 
thus becomes 



2M+1 



pisi\s)= n p{s^\st'\-^'-''-'^)pis^^\4-'\-^'-''-'^) 



with 



p( + l| + l,^)=p(-l 
p(-l| + l,A)=p( + l 

pi+l\-l,A)=p{-l 
p{-l\-l,A)=p{+l 
p{-l\ + l,B)=pi+l 
p{+l\ + l,B)=p{-l 

Introducing again the binary variables 



-l,i^) = l-7oAi 

+ 1, D) = 1 - -f2At 

- 1, B) = p{+l\ -l,C)= p{-l\ + 1,C)= 71 At 

- 1, B) = p(+l| + 1, C) = p{-l\ -1,C) = 1- 71 Ai 



<,2 = ^^, 6^,2 = ^^^ , 0<i^<M, 

we then find 

p{Si\s) = {i-por'"Po"a~pir'Pi'i^-P2r'PT 

where po = 7oAf,pi = jiAt,p2 = 72At, and mo,TOi,m2, no,ni,n2 are given by 

M 

mo = E {<'^2ibr'K + br%) + (1 - <)(! - a2)[(l - b^'K^ &i) + (1 - &r')(l - b^ 

i/=0 

Af 

no = 5] {<a^[&r'(i - &i) + &r'(i - &2)] + (1 - <)(i - «2)[(i - br'K + (i - for')fe2]} 

i'=0 
M 
771,1 

l/=0 
Af 
ni 

i/=0 
A/ 

m2 = 5] {(1 - <)(! - a'^W.-'b^i + hl-^ir^) + «i«2([(l - b\-^){\ b^) + (1 - 6r')(l - ^2)]} 

i/=0 

Af 



^1 = E {K(i - «2) + (1 - <)a2][fcr'&i + 6r'&2 + (1 - fcr')(i - ^d + (i - ^r')(i - ^d]} 

i/=0 

Af 

^1 = E {K(i - «2) + (1 - <)«2][&r'(i - &i) + br\^ &2) + (1 - &r')&i + (i - br'm} 



"■2 



i/=0 



(A17) 



(A18) 

(A19) 
(A20) 



E {(1 - <)(i - «2)K"'(i - &D + &2-'(i - ^2)] + <«2[(i - br'K + (1 - ^'r')fe2]} • (A21) 



with bi = af^ and 63 = a^^ (note that v is here an index and not an exponent). The first task is to prove Eq. (A3) 
as in the case A^ = 1. This amounts to proving that the ratio 



R: 



P{Si\S)PiS\Si) 



P{Si\S)PiS\Si) 
is equal to 1. After lengthy calculations we find that 

2 

J=0 



(A22) 



(A23) 
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with 



fco = ^2 
lo ^ h - 



fci 
2 
h 
2 ' 



(A24) 



M 



^1 = E {[^1^' - ^1 + ^2 + ' - ^2 + 2(6r'&r' - 6rfe2^)][«' + 44+' + (1 - <)(1 - <+') + (1 - «2)(1 - «2 + ')] 

+ K+1 - < + 4+' - "2 + 2K+'ar' - «l«2)][6lfer' + fc2&r' + (1 - &D(1 - &! + ') + (1 - &2)(1 - ^2 + ')]} 
M 

h - Eit^i^' - ^1 + ^2^' - ^2 + 2(6r'&2^+' - foi&2)]K(i - <+') + 4(1 - 4+') + (1 - <)<^' + (1 - 4)4+'] 






,1^+1 



^ _ aj + ar^ - 4 + 2(4+'4+' - «i4)][6i(i - 6r') + ^^2(1 - hT') + (1 - fei)&r' + (1 - &2)(i - ^r')]} 

(A25) 



with the convention a^ ^^ = ^12 ^-^d 6^ 2+' 



"^[(i 



Hence 



(1-P0)(1-P2)^ ^P0P2^ 



(A26) 



The two numbers fci and li are not zero and therefore R is not equal to 1 at this stage. However, one can easily 
convince oneself that fci and li are related to the number of switchings of the system during the time interval r 
(for instance, we see in the expression of fci that the first term inside brackets is zero if h'^'^^ — h\ and 63 +' = ^21 
which means that the state of the system has not changed in the corresponding time step Ai). Since this number of 
switchings remains finite in the continuous limit Ai — > 0, we then have 



n _ I- (1 - ll^t) iKi/.^|- fl -iti/-^ / /I X 

'-(1 -7oAt)(l -72Ai)-' '-7072 7o72 



fci/2r 7? ^ll/2 ^ I 7i_w,/2 



(A27) 



when At — > 0, and we conclude that the transition rates must satisfy Eq. (10) in order to have i? = 1. 




FIG. 7: (Color on line) Probabilities •p(A,t — T\A,t\A) (black curve) and p{A,—t;A,T — t\A) (red curve) for r = 200 when 
the relation between the rates, Eq. (10 1, is not satisfied. This case corresponds to the arbitrary choice 70 = 72 = 0.003965 and 



71 — 0.001516. One can see that time-symmetry does not hold in this case. 



The rest of the proof proceeds as in the case N = 1. To illustrate the above demonstration, we show in Fig. [7] an 
example where Eq. ( 10 1 is violated and time-symmetry does not hold. 
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Appendix B: Solution of t lie set of Eqs. (21) 



In this appendix we detail the solution of the set of linear differential equations ( 21 ). As in section III, all probabilities 
refer to the stationary state. 

The expressions of the matrices Ta and Tb are 



and 
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/ 



(Bl) 



(B2) 



where 7 = 70 + 72 — 271 . 

These matrices have 4 zero eigenvalues because the functions p{s',t — T;s,t|s°) are not linearly independent. This 
is a consequence of the conservation of probabilities, 



^p(s',i-r;s,f|s°) ^ 1 

s,s' 

^p(s',t-r;s,t|s")p(sO)=p(s) 
^p(s',i-r;s,t|s>(s°)=p(s'). 



(B3) 



(Of course, these relations may be used to reduce the size of the matrices.) The remaining 6 eigenvalues of Ta and 
Ffc, denoted ±Ayi^i,±AA,2, ±Ayi^3, and ±Xb,i,±Xb,2,^^b,3, respectively, are given by 



Xa,i = Y 7^ + 87? + 87072 

Aa,2 = 2 Wi^ + 167072 + 7] 

Aa,3 = 2 [^7^ + I67072 - 7] 



(B4) 



and 



Asa = y 7^ + 87? + 87072 
Ab,2= 2[V7' + 167?+7] 
Ab,3 = ^ [772 + 167? -7] 



(B5) 
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Therefore the two matrices have the same spectrum if the relation ( 10 ) 7072 
then simply denoted ±Ai,±A2,±A3 with 



jf is satisfied. The eigenvalues are 



Ai 



^1' + I67? 



A2 = ^{M + j) 
^3 = ^{>^i -7) 



and Eqs. (23) may be recast as 






(B6) 



(B7) 



where R^, Hb are diagonal matrices with (— Ai, Ai, — A2, A3, —A3, A2, 0, 0, 0, 0) and (— Ai, Ai, —A3, A2, — A2, A3, 0, 0, 0, 0) 
on the diagonal, respectively, and the matrices Ua,Ub read 
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p{D, r\A)^^ [-Ai + A2e^-^^ + X^e'^n 



(B8) 



(B9) 



(BIO) 



(Bll) 
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and 



p(B, r\B)^^ [Aie-^^^ + Aae-^^^ + A3] 



piC t\B) = -— [-Aie-^^^ + Aae-^i- + A3] (B12) 

with 

Ka = Ai(72 - 70) + [(70 + 72)A3 - 'ill]e-^'^ + [(70 + 72)A2 + i-/^]e^''' 

Kb = 2[A3 + 271 + (A2 - 27i)e-^i"] . (B13) 

As noted in section III, these nontrivial solutions of the self-consistent equations only exist when the rates satisfy Eq. 

diol. 



The stationary probability p{A) is most easily computed by setting t = r in the first equation (20 1 which expresses 
the condition of detailed balance between the states A and B in the stationary state 



This yields 



_ 1 p{A,r\B) 

^^'^^-2piA,r\B)+p{B,r\A)- ^^^'^ 



piA) . 1 — f^t^-^'^l ^ . (B15) 



In particular, the values of p{A) for t = and r — > 00 are 



and 



PoiA) ^ \ (B16) 

2(71 +72) 



.,^_ 1 47f + A2(7o + 72) , . 

^°°^^^ - 2 87? + A2(7o + 72) + 271 A3 ■ ^ '> 
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